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PART  I 


8ENDING  OF  THICK  BEAMS  OF  BIMODULUS  MATERIALS 
A.D.  Tran 

Exxon  Production  Research 
Houston,  Texas 

and 

C.W.  Bert 

University  of  Oklahoma 
Norman,  Oklahoma 

Abstract  -( The  1  iterature  on  bending  of  beams  made  of  bimodulus  materials 
(which  have  one  value  of  the  elastic  modulus  in  tension  and  another  in  com¬ 
pression)  is  limited.  All  of  the  works  known  to  the  present  investigators 
are  restricted  to  beams  with  natural  boundary  conditions  and  subjected  to 
specific  distributions  of  normal  pressure  only.  In  the  present  work,  the 
transfer-matrix  approach  is  used  to  determine  the  small -deflection  static 
behavior  of  bimodulus  beams,  including  transverse  shear  deformation.  The 
neutral  surface,  i.e.,  the  locus  of  points  having  zero  axial  normal  strain, 
may  vary  linearly  within  each  element.  The  effects  of  axial  load  and  non¬ 
natural  boundary  conditions  are  considered.  As  a  basis  for  comparative  evalua¬ 
tion,  exact  closed-form  solutions  are  also  presented  for  special  cases  in  which 
the  neutral -surface  location  is  constant  along  the  beam  axis.  Results  are 
compared  between  the  two  solution  methods  and  are  found  to  give  good  agreement. 

INTRODUCTION 

As  early  as  1864,  St.  Ve.nant  [1]  recognized  that  certain  actual  materials 
have  different  elastic  behavior  when  they  are  loaded  in  tension  as  compared 
to  compression.  See  also  [2-4].  However,  the  concept  of  a  bimodulus  material, 


i.e.,  a  bilinear  material  having  different  moduli  in  tension  and  in  compression 
was  not  originated  until  1941  hy  Timoshenko  [3],  who  considered  the  flexural 
stresses  in  such  a  material  undergoing  pure  bending.  The  effective  modulus 
for  stiffness  of  such  a  beam  in  pure  bending  was  given  by  Marin  [3].  The 
bimodulus  concept  was  rediscovered  and  extended  to  two-dimensional  materials 
by  Ambartsumyan  [6 ]  in  1965.  Subsequently,  there  have  been  many  analyses  of 
two-dimensional  bimodulus  material,  but  this  topic  is  beyond  the  scope  of 
the  present  work. 

Work  analogous  to  that  of  Timoshenko  [3]  or  Marin  [5]  on  static,  small  - 
deflection  bending  of  Bernoull i -Euler  beams  of  homogeneous,  bi modulus  mate¬ 
rial  was  presented  in  [7-15].  A  so-called  no-tension  material,  which  has  no 
resistance  to  tension  but  a  finite  elastic  resistance  to  compression,  is  a 
special  case  of  bimodulus  material  applicable  to  brittle  materials,  such  as 
many  ceramic  materials.  Buckling  of  columns  of  no-tension  material  was  ana¬ 
lyzed  in  [16-18].  Column  buckling  of  general  bimodulus  material  was  con¬ 
sidered  in  [19,20].  Large  static  deflections  of  beams  of  bimodulus  material 
were  treated  in  [21,22].  Transverse  shear  effects  on  bimodulus  beams  were 
first  treated  by  Kamiya  [23].  Bimodulus  action  in  viscoelastic  beams  was 
considered  by  Nachlinger  and  Leininger  [24]. 

Special  studies  applicable  to  bending  of  specific  kinds  of  nonhomogeneous- 
material  beams  were  carried  out  by  some  investigators.  For  example,  Vierling 
and  Scheele  [25]  were  concerned  with  cord-rubber,  Bilida  [26]  and  Phan-Thien 
[27]  with  fiber  composites,  Carlsson  et  al .  [28]  with  paperboard,  Swift  [29] 
with  reinforced  concrete,  and  [30-33]  with  laminates. 

The  present  work  is  believed  to  be  the  first,  in  the  context  of  bimodulus 
beams,  in  the  following  respects: 
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1.  Consideration  of  complicated  and  concentrated  loadings,  the 
latter  causing  the  neutral -surface  location  (and  thus  the  beam 
stiffnesses)  to  change  along  the  length  of  the  beam. 

2.  Application  of  the  transfer-matrix  method. 

GOVERNING  EQUATIONS  FOR  SMALL  DEFLECTIONS  OF  A  BIMODULUS 
BEAM  ACCOUNTING  FOR  TRANSVERSE  SHEAR 


Consider  a  beam  of  thickness  h  and  length  i.  The  origin  of  a  Cartesian 
coordinate  system  is  located  on  the  central  axis  with  the  z-axis  being  normal 
to  the  central  axis  (see  Fig.  1).  The  stress  resultants  and  stress  moment, 
each  per  unit  width,  are  defined  as: 


rh/2 

rh/2 

]  (V  Txz)dz 

;  M  = 

z,xdz 

(1) 

-h/2 

-h/2 

where  ax  and  txz  are  respectively  the  axial  normal  stress  and  the  transverse 
shear  stress.  The  theory  developed  by  Yang,  Norris,  and  Stavsky  [34]  assumes 
the  following  displacement  field: 

U(x,z)  =  u(x)  +  z^(x)  ;  W(x,z)  =  w(x)  (2) 


where  U  and  W  are  the  displacement  components  in  the  x  and  z  directions, 
respectively,  u  and  w  are  the  corresponding  displacements  at  the  midplane, 
and  -I/  is  the  bending  slope, 

The  constitutive  relation  for  a  bimodulus  beam  with  provision  for  a 
shift  of  the  neutral  surface  (due  to  different  properties  in  tension  and 
compression )  can  be  written  as  follows  [35]: 
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Here,  differentiation  is  denoted  by  a  comma,  i.e.,  (  ),  s  d(  )/dx.  The 
symbols  A,  B,  D,Sdenote  the  respective  extensional,  flexural-extensional 
coupling,  flexural,  and  thickness  shear  stiffnesses  defined  by 


(A,B,D) 


h/2 

Q^O ,z,z2)dz 

-h/2 


-h/2 


(k  =  t,c) 


(4) 


The  quantity  K2  is  a  shear  correction  coefficient  which  is  generally  taken 
to  be  5/6  for  static  loading  of  rectangular-section  beams;  is  E^  for  a 
compact-section  beam  or  E^kV(l  -  vi2^v21^)  f°r  a  wide,  thin  beam; 
is  the  longitudinal -transverse  shear  modulus;  and  k=t  for  tensile-strain 
regions  or  k  =  c  for  compressive-strain  regions. 

One  can  write  the  equations  of  equilibrium  as  follows: 


N,x  *  0  i  <?.x  =  ‘  ;  M,x  -  Q  =  0  (5) 

Here,  q(x)  is  the  transverse  distributed  loading. 

Substituting  eqns  (3)  into  eqns  (5),  we  obtain  the  following  equations 
of  equilibrium  in  terms  of  generalized  displacements 

(Au,x  +  B^,x),x  =  0 

[S(w,x  +  *)],x  =  -  q(x)  (6) 

(Bu,x  +  Dli; ,x) , x  -  S(w,x  +  tp)  a  0 

The  expression  for  neutral -surface  location  [36],  zn>  can  be  written  as: 

c  =  u ,  +  z  ii,  =  0 
x  ’x  n^’x 

or 

zn  ‘  -  “•x/*,x  <7> 


ANALYSIS:  CLASSICAL  APPROACH 


For  cases  when  the  neutral-surface  locations  are  constant  along  the 
beam  axis  (x-direction) ,  eqns  (6)  become 

Au’xx  +  B*’xx  =  0 

A(w,xx  +  *.x)  =  -q(x)  (8) 

Bu’xx  +  D*»xx  ‘  S(w’x  +  =  0 

Eqns  (8)  constitute  a  sixth-order  system  of  linear  ordinary  differential 
equations  with  constant  coefficients.  The  general  solution  can  be  written 
as  [37]: 

3B 

u(x)  ■  dx  +  d2x  +  -J-  c4x2  +  Up 
2 

*(x)  =  -c2  +  6^B--S'a~^  C4  '  2c3X  -  3c4x2  +  *p  (9) 

w(x)  s  Cl  +  C2X  +  C3x2  +  c4X3  +  Wp 

where  Cj,  c2»  c3,  c4  are  arbitrary  constants  and  u  ,  ■>  ,  w  are  particular 

r  r  r 

solutions.  The  arbitrary  constants  are  determined  by  the  boundary  conditions 
of  the  beam  at  x  =  0  and  x  =  L. 

The  clamped  boundary  conditions  at  x  =  0  and  x  =  l  are  given  as 

u(0)  =  uU)  =  0  ;  j(0)  -  -p(£)  =  0  ;  w(0)  =  w(£)  =  0  (10) 

The  clamped  boundary  conditions  at  x  =  0  and  free  boundary  conditions 

at  x  *  i  are  given  as 

u(0)  *  N (0  »  0  ;  *(0)  -  M(e)  -  0  ;  w(0)  =  Q(i)  =  0  (11) 

The  values  of  constants  clt  c2,  c3,  c4,  d2,  d2  for  the  clamped-clamped  and 

clamped-free  boundary  conditions  are  as  follows: 
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sq0 

up{x)  s  ^AB-~b3)  cos  aX 

Aq0 

Vx)  =  '  IWrW  cos  aX 


(15) 


q  i  ^ 

WP(x)  =  $  [S  +  7^AD“B~)]  S':‘l  aX 


For  normal  load  q(x)  =  qQcos  ax,  the  particular  solutions  are 

Bq 

up(x)  =  “  ^(AD  -  B^)  S1n  aX 
Aq 

*pM  =  a^(AD  -°B l)  Sln  aX 


06) 


wp(x)  =  rr  Ct  + 


$  [S  +  a^(AD  -  B") 3  C0s  aX 


ANALYSIS:  TRANSFER-MATRIX  APPROACH 


A  common  type  of  system  occurring  in  engineering  practice  consists  of 
a  number  of  elements  linked  together  end  to  end  in  the  form  of  a  chain. 
Well-known  examples  are  continuous  beams,  turbine-generator  shafts,  etc. 

The  transfer-matrix  approach  is  ideally  suited  to  such  systems  because  only 
successive  matrix  multiplications  are  necessary  to  fit  the  elements  together. 

In  the  transfer-matrix  techniques  [38,39],  the  beam  is  divided  into  Ng 
elements  of  mass  m  and  each  mass  m  is  assumed  to  be  concentrated  at  the  mass 
center  of  the  element.  The  locations  of  these  point  masses  are  referred  to 
as  stations.  Consecutive  stations  are  separated  by  massless  fields  which 
contain  all  of  the  stiffnesses  of  the  system.  Consequently,  the  model  of 
the  beam  shown  in  Fig.  2  consists  of  two  half  fields,  aE/2,  one  at  each  end 
of  the  beam,  and  N  stations  separated  by  N$  -  1  full  fields  of  length, 
where  aE  =  E/Ns  and  i  is  the  total  length  of  the  beam.  The  transfer  matrix 


* 


til 


transfers  the  displacements  (u,w,'Jj)  and  the  forces  (M,Q,N)  from  the  left 
side  of  the  field  or  station  to  the  right  side  of  the  same  field  or  station. 

Figure  3  shows  the  flexural  configuration  of  the  i^  station.  From 
continuity  and  static  equilibrium  considerations  of  Fig.  3,  the  following 
conditions  are  derived: 


u1: 

mr 

=  ml 
1 

R 
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qt 
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The  above  equations  can  be  written  in  matrix  notation  as  follows: 


Matrix  eqn  (18)  can  be  written  in  more  compact  notation  as 

[SlJ  =  [TglflSlf  09) 

where  the  definitions  of  each  term  can  be  ascertained  easily  by  direct 
comparison  between  eqns  (18)  and  (19).  The  matrix  [T  ]  is  known  as  the 
station  matrix. 

Figure  4  illustrates  the  flexural  configuration  for  the  itn  massless 
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beam  element  (field).  Here,  [S]^  (station  i)  =  [ S ] V  (field  i)  and  [S]^ 
(field  i)  =  [S]l:+1  (station  i+1).  From  equilibrium  considerations  the 
following  expressions  are  evident: 

NW  ■  «?  !  8^  =  5?  -  K,  ;  <t,  -  -  K, 


where 


rtJL  fiu. 

Kq  =  q(c)d$  ;  Km  =  (  SqU)<U 
0 


r 

0 


Eqns  (3)  can  be  inverted  to  become: 


(20) 


(21) 


u,. 
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- 

(22) 


w,x  =  (Q/S)  -  * 


From  eqn  (22),  one  can  write  the  following  equations  by  taking  the  averages 
of  the  respective  moments  and  axial  forces  at  both  ends  of  the  field: 

B(M^+1  +  M*)  0(N|[+1  +  N?) 


"HI  ’  “?  +  f  '■  2 


-] 


(23) 


,1  _  ,R  bJL  , 

*1+1  '  *i  -  T  [ 


A(M*T+^  +  M?)  B(N1r+1  +  N?) 


where  ■<  =  B2  -  AD. 

Making  use  of  eqns  (20),  one  can  write  equations  (23)  as: 

V!  •  uj  +  (§^)M*  +  [BMII.jqR  .  (M)n*  .  (M)Km 
♦i+1  =  *?  -  (~)M?  -  +  (M)nR  +  (M)Km 


(24) 


Similarly,  eqns  (22),  (20),  and  (23),  can  be  combined  to  yield  the  following: 


Values  for  K  and  K  for  several  different  loading  conditions  are  listed  in 
Table  1.  q  171 


(25) 


which  can  also  be  written  as 

[S]-+1  =  CTj] CS] i  (27) 

The  matrix  [T.J  is  known  as  the  field  matrix. 

The  station  and  field  matrices,  eqns  (19)  and  (27),  can  be  used  to  give 
the  relationship  of  the  state  vectors  at  the  two  ends  of  the  beam  by  system¬ 
atically  applying  the  matrices  to  each  station  and  field  as  follows: 


where  M/2  is  the  length  of  each  of  the  half  fields  at  the  ends  of  the  beam 
and  M  is  the  length  of  each  of  the  whole  fields.  Eqn  (28)  can  be  written 
in  compact  form  as 


[S]N^  »  (R][S]q  (29) 

Given  the  boundary  conditions,  eqn  (29)  can  be  solved  for  the  state  vectors 
at  the  ends,  which  are  [S]  and  [S]N+1;  see  Appendix  C  of  [40]. 

Finally,  the  rest  of  the  state  vectors  are  calculated  as  follows: 


o 


where  the  superscripts  L  and  R  refer  to  the  left  and  right  sides  of  the 
station,  respectively. 

Values  of  beam  stiffnesses  A,  B,  D,  S  are  needed  in  order  to  solve  the 
governing  equations.  However,  for  bimodulus  beams,  neutral-surface  positions 
are  required  for  the  computation  of  A,  B,  D,  S.  Study  of  the  expression  for 
neutral -surface  position  in  eqn  (7)  indicates  that,  except  for  a  few  special 
cases  when  neutral -surface  positions  are  known  constants,  the  governing 
equations  must  be  solved  before  we  can  compute  neutral-surface  positions. 

Our  approach  is  first  to  assume  a  set  of  neutral-surface  positions 
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along  the  beam  length,  then  solve  the  governing  equations  for  displacements 
u,  w,  i>,  and  the  derivatives  u,  and  ■  With  known  values  of  u,  and  -j >,  , 

XX  XX 

compute  neutral -surface  positions  using  eqn  (7).  If  the  assumed  and  computed 
set  of  neutral -surface  positions  are  in  close  agreement,  the  problem  is 
solved.  Otherwise,  assume  a  different  set  of  neutral -surface  positions  and 
repeat  the  calculations. 


NUMERICAL  RESULTS 

Numerical  results  were  obtained  to  verify  the  validity  of  application  of 
the  transfer-matrix  technique  to  bimodulus-beam  problems  and  to  study  the 
behavior  of  beams  made  of  aramid-cord  rubber  (a  bimodulus  material).  Unless 
stated  to  the  contrary,  the  material  properties  and  beam  dimensions  listed 
in  Table  2  are  employed  in  all  of  the  problems  considered,  which  are  summarized 
in  Table  3.  Values  of  the  constants  used  in  Table  3  are  as  follows: 

Mq  =  =  0.1  lb-in.  ;  =  1.0  lb  ;  qQ  =  q,  =  0.1  lb/in. 

The  transfer-matrix  model  of  the  beam  includes  twenty-five  elements, 
each  of  length  0.32  in.  A  value  of  5/6  was  used  for  the  shear  correction, 
coefficient,  K2.  All  of  the  computations  were  carried  out  on  an  IBM  370 
Model  158J  computer. 

Response  of  an  aramid-cord  rubber  beam  subjected  to  various  loadings  is 
investigated  in  Cases  1  through  8.  Bending  moment  M,  shear  force  Q,  axial 
force  N,  and  the  corresponding  displacements  including  axial  elongation 
(contraction)  u,  transverse  deflection  w,  bending  slope  v,  together  with 
neutral -surface  positions  zn>  are  computed.  Results  are  presented  in  Tables 
4  through  11  for  Cases  1  through  8,  respectively. 
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The  boundary  conditions  of  Cases  1  through  4  are  fairly  simple.  We 
have  been  able  to  develop  both  closed-form  and  transfer-matrix  solutions 
for  these  cases.  The  good  agreement  between  closed- form  and  transfer-matrix 
results  can  be  seen  in  Tables  5  through  8.  For  the  majority  of  the  results, 
the  error  is  less  than  2%,  although  for  a  few  result  values  which  are  close 
to  zero,  the  percentage  error  can  go  up  to  5%.  Actually,  accuracy  of  the 
transfer-matrix  results  can  be  improved  by  increasing  the  number  of  elements 
into  which  the  beam  is  divided.  The  number  of  elements  used  in  this  study 
is  twenty-five  which  is  judged  sufficient  to  give  reasonably  good  results. 
Support  of  the  judgment  is  shown  in  Fig.  5. 

Behavior  of  cl amped-free  bimodulus  beams  subjected  to  uniformly  distri¬ 
buted  and  sinusoidally  distributed  loads  is  investigated  in  Cases  1  and  2. 
Behavior  of  bimodulus  beams  which  are  also  subjected  to  uniformly  distributed 
and  sinusoidally  distributed  loads  but  with  clamped-cl amped  boundaries  is 
investigated  in  Cases  3  and  4.  The  case  of  simply-supported  boundaries  and 
sinusoidal  load  can  also  be  solved  in  closed  form  for  a  beam.  However,  this 
condition  is  not  discussed  here  because  it  is  a  special  one-dimensional  case 
of  the  plate  problem  investigated  in  [36].  For  all  four  Cases  1  through  4, 
the  axial  force  N  is  zero.  Neutral -surface  positions  zn  are  either  constant 
throughout  or  piecewise  constant.  Values  of  zn  for  Cases  1,  2,  and  3  are 
-0.2672,  0.2672,  and  0.2672  in.,  respectively.  The  value  of  zn  for  Case  4 
is  0.2672  in.  when  0<_x<  i/2  and  is  -0.2672  in.  when  z/2  <  x  <_  i.  Actually, 
zero  axial  force  always  leads  to  constant  or  piecewise  constant  neutral- 
surface  position;  this  can  be  shown  using  eqn  (5.3)  from  Appendix  B  and 
setting  N  =  0  and  M  t  0: 

z„  =  B/A  =  constant  (31 ) 

n 


Constant  or  piecewise-constant  neutral -surface  position  greatly  reduces  the 
difficulty  in  development  of  closed-form  solutions  for  eqns  (6),  because 
eqns  (6)  can  then  be  reduced  to  a  much  simpler  form,  eqns  (8). 

If  a  practicing  engineer  does  not  have  any  access  to  the  technique  of 
analyzing  bimodulus  structural  elements  as  presented  here,  the  simplest 
thing  he  would  do  is  to  treat  the  bimodulus  beam  as  an  ordinary  beam  with 
Young's  modulus  taken  to  be  the  arithmetic  average  of  the  tensile  and  com¬ 
pressive  Young's  moduli  of  the  bimodulus  material.  The  question  of  how  good 
is  this  approach  then  arises;  the  answer  is  presented  in  Figs.  6  and  7.  The 
distribution  of  maximum  tensile  and  maximum  compressive  normal  stresses  at 
the  outermost  layer  of  the  beam  subjected  to  the  condition  of  Cases  1  and  2 
are  presented  for  both  approaches  in  Figs.  6  and  7.  In  a  similar  fashion, 
deflections  are  compared  in  Fig.  8.  As  shown  in  Figs.  6  and  7,  the  average- 
modulus  approach  results  in  a  maximum  normal  stress  which  is  approximately  10% 
of  the  actual  maximum  stress  (resulting  from  the  bimodulus  approach)  for  both 
Cases  1  and  2.  Figure  8  shows  that  a  maximum  deflection  approximately  10%  of 
the  actual  deflection  (based  on  the  bimodulus  approach)  also  results  from  the 
average-modulus  approach.  In  short,  design  based  on  the  average-modulus 
approach  would  likely  lead  to  premature  failure  of  the  structural  elements. 

The  technique  developed  in  this  work  therefore  would  be  valuable  to  design 
engineers. 

To  compare  the  results  of  both  the  transfer-matrix  analysis  and  finite- 
element  analysis  with  the  closed-form  solution,  Figs.  9  and  10  are  presented 
for  Case  3.  The  finite-element  results  were  provided  by  Dr.  J.N.  Reddy  [41] 
and  were  obtained  as  a  special  case  of  a  bimodulus  plate  [42].  There  were 
twenty-two  elements  along  the  length  of  the  beam.  Even  for  the  exaggerated 


scale  used  in  these  figures,  it  was  not  possible  to  show  a  difference  between 
the  transfer-matri x  (TM)  and  closed-form  (CF)  solutions.  The  finite-element 
values  were  generally  smaller  than  the  CF  and  TM  values,  except  for  bending 
moment,  for  which  all  three  methods  agreed  very  closely.  It  was  not  possible 
to  compare  the  computation  times  involved,  since  the  results  were  run  on 
different  machines.  However,  in  the  past,  it  has  been  the  general  experience 
that  the  TM  technique  is  more  efficient  [43]. 

Since  closed-form  solutions  are  not  available  for  the  complicated 
boundary  conditions  considered  in  Cases  5  through  8,  only  transfer-matrix 
results  are  presented  for  these  cases.  The  effect  of  concentrated  moment, 
axial  force,  and  transverse  force  applied  at  the  free  boundary  on  static 
behavior  of  a  clamped-free  bimodulus  beam  are  investigated  for  two  distributions 
of  loading  (uniform  and  sinusoidal)  in  Cases  5  and  6.  The  effect  of  concentra¬ 
ted  moments  applied  at  both  boundaries  on  static  behavior  of  a  bimodulus 
beam  subjected  to  the  boundary  condition  of  hinged-hinged  with  axial  constraint 
are  investigated  for  uniformly  distributed  and  sinusoidal  loads  in  Cases  7  and 
8.  As  shown  in  Tables  8  through  11,  the  axial  forces  N  induced  by  the  con¬ 
centrated  moment  and  forces  at  the  boundaries  are  constant  but  not  zero;  and 
neutral -surface  positions  z^  vary  drastically  along  the  beam  length.  In  fact 
a  look  at  eqn  (8.3)  tells  us  that  if  axial  force  N  is  nonzero,  one  would  expect 
zn  generally  not  to  be  constant.  Although  for  beams  made  of  ordinary  (not 
bimodulus)  material,  the  neutral  surface  always  coincides  with  the  geometric 
midplane,  the  neutral  surface  for  a  bimodulus  beam  is  not  only  far  different 
from  the  midplane  but  also  assumes  various  shapes  depending  on  the  conditions 
of  boundaries  and  leadings.  Shapes  of  neutral  surfaces  for  Cases  5  through  8 
are  illustrated  in  Fig.  11  through  14. 
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CONCLUDING  REMARKS 

In  this  study,  the  transfer-matrix  model  of  a  bimodulus  beam  was 
developed  based  on  shear-deformable-beam  theory.  Transfer-matrix  results 
were  compared  with  and  were  in  close  agreement  with  closed-form  solutions 
developed  herein  for  all  four  combinations  of  uniformly  distributed  and 
sinusoidally  distributed  load  with  clamped-clamped  and  clamped-free 
boundary  conditions. 

Although  closed-form  solutions  are  available  only  for  a  number  of 
simple  boundary  conditions,  the  transfer-matrix  model  offers  a  solution 
method  for  a  wide  variety  of  boundary  and  loading  conditions.  Utilizing 
the  model,  behavior  of  a  bimodulus  beam  was  investigated  by  computing 
and  observing  deflections  and  neutral -surface  positions  under  fairly  com¬ 
plicated  conditions  of  loading  and  boundary.  It  was  found  that  concentrated 
forces  and  moments  applied  at  the  boundaries  have  significant  influence  on 
neutral -surface  positions;  more  specifically  they  can  cause  the  neutral- 
surface  positions  (thus  stiffnesses)  to  change  drastically  along  the  beam 
length.  As  illustrations,  plots  of  neutral -surface  positions  were  presented 
for  some  of  the  investigated  cases.  Numerical  results  of  displacements  and 
stress  distributions  were  al so  presented. 

A  distinguishing  feature  of  the  model  developed  in  this  work  is  numer¬ 
ical  efficiency.  It  requires  much  less  computer  storage  than  the  variational 
finite-element  method.  Also,  in  a  comparison  to  an  exact  closed-form  solu¬ 
tion  presented  here,  the  transfer-matrix  results  were  more  accurate  than  the 
variational  finite-element  results. 
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APPENDIX  A.  EFFECT  OF  NEUTRAL-SURFACE  LOCATION  ON  THE  BEAM  STIFFNESSES 

It  is  mathematically  convenient  in  the  solution  of  eqns  (6)  to  reduce 

these  equations  to  an  algebraic  form  not  including  the  integral  forms  of  the 

A,  B,  D,  and  S.  Therefore,  the  expansion  of  these  terms  into  their  form  in 

k  k 

terms  of  the  two  Young's  moduli  E  and  two  shear  rigidities  G  is  completed  in 
the  following  manner. 

Refer  to  Fig.  A.l  for  a  typical  stress  distribution  for  the  case  of  the 
neutral -surface  location  greater  than  zero;  where  z  ,  the  neutral -surface 
location,  is  measured  positively  downward  from  the  geometric  midplane.  The 
top  of  the  beam  is  considered  to  be  -h/2  and  the  bottom  of  the  beam  is  h/2. 

The  beam  stiffnesses  are  defined  as  follows: 


fh/2  V 

E  az  , 

-h/2 

rh/2  k 

B  =  E*z  dz 

■'-h/2 

(A.l) 

fh/2  k  2 

E  zaz 

-h/2 

■  S!J 

rh/2  . 

G*dz 

-h/2 
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1/  1/ 

where  E  and  G  are  the  Young's  moduli  and  shear  moduli  respectively . 
Superscript  k  can  be  either  t  or  c,  which  represent  tension  or  compression 
properties  respectively. 

Convex  Downward  Bending 

In  convex  downward  bending,  the  top  layer  of  a  beam  is  in  compression 
and  the  bottom  layer  in  tension. 

The  first  of  eqns  (A.l)  can  be  expanded  in  the  following  form: 


,zn  h/2 

A  =  f  Ecdz  +  f  E*dz  =  [z„  +  (h/2) ]EC  +  [(h/2)  -  z^E* 


-h/2 


n 


n 


or 


A  =  (h/2) (EC  +  E*)  +  (EC  -  Et)zr 


(A. 2) 


In  similar  form,  the  second  of  eqns  (A.l)  can  be  expanded  to: 


f  n 

3  =  Ecz  dz  + 

-h/2 


zn 


or 


B  =  -  (h2/3) (Ec  -  Et)  +  (1/2) (Ec  -  E^z* 


(A. 3) 


Also 


zn 

f  n  r  2 
D  =  i  Ecz  dz  + 

-h/2 


h/2  3  3 

Etz2dz  =  (1/3) (z2  +  ^-)EC  +  (1/3)(V  -  0Et 


8  n 


or 


D  =  (h3/24)(Ec  +  E*)  +  (1/3) (EC  -  Et)z^ 


(A. 4) 


h/2 


S  = 


Gcdz  + 


Gtdz  =  [zn  +  (h/2) ]GC  +  [(h/2)  -  zn]Gt 


-h/2 


and 
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S  =  (h/2) (GC  +  G1)  +  (GC  -  Gt)zn 
Concave  Downward  Bending 

In  concave  downward  bending,  the  top  layer  of  beam  is  in  tension  and 
the  bottom  layer  in  compression. 

The  first  of  eqns  (A.l)  can  be  expanded  in  the  following  form: 


Etdz  +  f  ECdz  =  Et[zn  +  (h/2)]  +  Ec[ (h/2)  -  zp] 


A  =  (h/2) (EC  +  E1)  -  zn(EC  -  Et) 


With  a  similar  approach,  the  second  of  eqns  (A.l)  can  be  expanded  to: 
h/2  *  ,  - 

B  -  {  Etz dz  +  J  Eczdz  =  V  (zn  "  T}  +  T  (T  '  zn> 

-h/2  zn 


B  =  (h2/8)(EC  -  El)  -  (1/2) (Ec  -  E^z2 


0-("  E'zdz*!  E/dz  =  T(zn  +  T>  +T(T-  zn> 

-h/2  zn 


D  -  (h3/24)(EC  +  Efc)  -  (1/3) (EC  -  E^z2 


Gtdz  +  f  Gcdz  =  Gt(zn  +  |)  +  Gc(|  -  zR) 


S  =  (h/2) (GC  +  Gt)  -  (GC  -  Gt)z 


23 


APPENDIX  8.  METHOD  OF  COMPUTING  NEUTRAL -SURFACE  LOCATION  FROM  MOMENT 
DIAGRAM  AND  AXIAL  FORCE 

Under  certain  combinations  of  loading  and  boundary  conditions,  the 
neutral -surface  location  varies  drastically  along  the  beam  length  (x-axis). 
Closed-form  solutions  currently  exist  only  for  cases  in  which  the  neutral- 
surface  location  is  constant  along  the  beam  length.  However,  the  transfer- 
matrix  technique  can  be  utilized  to  analyze  cases  of  varying  neutral -surface 
location.  Evidently  some  benchmark  is  needed  to  verify  the  latter  technique's 
application  to  bimodulus  problems.  It  is  a  fairly  simple  task  to  derive  the 
moment  diagram  from  boundary  conditions  and  applied  axial  force.  Results  of  this 
method  have  then  been  checked  against  those  of  transfer-matrix  analysis; 
good  agreement  was  obtained. 

The  constitutive  relation  for  a  bimodulus  beam  with  provision  for  a  shift 
of  neutral -surface  location  (due  to  different  properties  in  tension  and  com¬ 
pression)  has  been  written  in  eqn  (22)  as  follows: 


B2 -  AD 


-D  B 
B  -A 


From  the  kinematics  of  deformation,  the  neutral-surface  location  has 
been  derived  and  expressed  in  eqn  (7)  as  follows: 

zn  =  '  u’>’x  (E 

Replacing  u°x  and  by  the  expressions  in  eqn  (B.l),  one  can  rewrite  eqn 
(8.2)  as 


Zn  =  (8M  -  DN ) / ( AM  -  BN) 


(8.3) 
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Study  of  eqns  (A. 2)  -  (A. 4)  and  (A.S)  -  (A. 8)  indicates  that  the 
expressions  for  stiffnesses  A,  3,  D  can  be  written  as 

A  '  al  *  a22n  ;  8  =  b2  *  b22i  ;  “  *  dl  +  d:zn  <8'4> 

where  for  concave  bending 

a:  =  (h/2) (Ec  +  Efc)  ;  a2  =  Ec  -  E1* 

bL  =  -(h2/8)(Ec  -  E*)  ;  b2  =  (1/2)(EC  -  Et)  (B.S.a) 

dt  -  (h3/24)(Ec  +  Et)  ;  d2  =  (1/3) (EC  -  E1) 

and  for  convex  bending 

ax  =  (h/2) (Ec  +  Et)  ;  a2  «  -  (Ec  -  EC) 

bL  =  (h2/8)(Ec  -  Et)  ;  b2  »  -  0/2) (Ec  -  Et)  (B.S.b) 

d,  =  (h3/24)(EC  +  Efc)  ;  d2  =  -  (1/3)(EC  -  E4) 

Replacing  stiffnesses  A,  8,  0  by  the  expression  in  eqns  (B.4),  one  can  write 
eqn  (B.3)  as 

M'2n2<a2-b2>  ♦Vn-'V  “  ‘  “iVt1"  (B'6) 

Once  moment  M  and  axial  force  N  are  known,  eqn  (B.5)  can  be  solved  for  zn . 
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Fig.  1 


Cartesian  coordinate  system  for  beam. 


elds  a£,  and  two  half  fields 


Fig.  4  Flexural  configuration  of  ith  field. 
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ordinary  beam  having  Young's  modulus  which  is  the  arithmetic 
average  of  tensile  and  compressive  Young's  moduli  of  aramid- 
cord  rubber  (Case  2). 
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E^l  REGION  IN  COMPRESSION 
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Fig.  11.  Neutral -surface  shape  of  an  aramid-cord 
rubber  beam  (Case  5). 


Fig.  12.  Neutral -surface  shape  of  an  aramid-cord 
rubber  beam  (Case  6). 


Fig.  13'.  Neutral-surface  shape  of  an  aramid-cord 
rubber  beam  (Case  7). 
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Table  5.  Comparison  between  closed-form  and  transfer-matrix  solutions 
for  an  aramid-cord  rubber  beam  (Case  2) 


denotes  closed-form  solutions;  TM  signifies  transfer-matrix  solution. 
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CF  denotes  closed-form  solutions;  TM  signifies  transfer-matrix  solution 


Table  8.  Forces,  displacements,  and  neutral -surface  positions 
of  an  aramid-cord  rubber  beam  (Case  5),  as  analyzed 
using  the  transfer-matrix  method* 
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orces,  displacements,  and  neutral -surface  positions 
f  an  aramid-cord  rubber  beam  (Case  6),  as  analyzed 
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Table  10.  Forces,  displacements,  and  neutral -surface  positions 
of  an  aramid-cord  rubber  beam  (Case  7),  as  analyzed 
using  the  transfer-matrix  method 
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NOTE:  Bending-stretching  coupling  caused  by  bimodulus  action  induces  an 
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Bending-stretching  coupling  caused  by  bimodulus  action  induces  an 
axial  force  N  =  -0.4004  It)  (compressive). 
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PART  II 

TRANSIENT  RESPONSE  OF  A  THICK  BEAM 
OF  BIMODULAR  MATERIAL 

C.W.  Bert 

School  of  Aerospace,  Mechanical  and  Nuclear  Engineering 
The  University  of  Oklahoma 
Norman,  Oklahoma,  U.S.A. 

and 

A.O.  Tran 
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Houston,  Texas,  U.S.A. 

SUMMARY 

Certain  materials  have  different  elastic  behavior  when  they  are  loaded  in 
tension  as  compared  to  compression.  As  an  engineering  approximation,  they 
are  usually  modeled  as  a  bimodular  material,  i.e.,  a  bilinear  material 
having  different  Young's  moduli  in  tension  and  in  compression.  All  of 
the  previous  analyses  of  bimodular  beams  known  to  the  present  investigators 
have  been  concerned  with  either  static  loading  or  harmonic  vibration.  Thus, 
the  present  work  is  believed  to  be  the  first  to  consider  transient  response 
of  such  beams.  The  transfer-matrix  method  is  used  to  discretize  spatially, 
while  the  timewise  discretization  is  accomplished  by  use  of  the  Newmark 


beta  method. 
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INTRODUCTION 

There  is  an  extensive  body  of  literature  on  the  static  behavior  of 
bimodular  beams,  going  back  as  far  as  Timoshenko  in  1941. 1  These  works 
were  recently  discussed  in  connection  with  a  transfer-matrix  static  analysis 
of  bimodular  beams  by  the  present  investigators.  Examples  of  material 
having  such  behavior  are  cord-rubber,  soft  biological  tissues,  paperboard, 
and  reinforced  concrete. 

The  only  dynamic  analyses  of  bimodular  beams  known  to  the  present 
investigators  are  the  works  of  Khachatryan3  in  1967,  Galoyan  and  Khachatryan4 
in  1978,  and  Tran5  in  1981.  All  of  these  involved  either  free  or  sinuso¬ 
idally  forced  vibrations. 

The  present  work  uses  the  transfer-matrix  method0,  which  was  shown  in 
Ref.  2  to  be  both  more  accurate  and  computationally  more  efficient  than  the 
finite-element  method  for  a  staticly  loaded  bimodular  beam.  The  neutral  - 
surface  position,  which  is  the  boundary  between  the  tensile  and  compres¬ 
sive  regions  of  a  cross  section  is  permitted  to  vary  in  a  piecewise  linear 
fashion  along  the  beam  length.  The  beam  is  modeled  as  a  Timoshenko  beam, 
i.e.;  both  transverse  shear  deformation  and  rotatory  inertia  are  considered. 

Numerous  approaches  have  been  used  to  handle  transient  response  prob¬ 
lems;  these  have  included  various  transform  techniques',  the  so-called 
direct-analysis  method9,  the  modal  method9,  and  various  numerical-integration 
schemes.  Among  the  latter  are  the  central -difference  method10,  the  Houbolt 
method11,  the  Newmark  beta  method1^,  the  Wilson  theta  method13,  and  stiffly 
stable  methods.14  For  simplicity  and  efficiency,  the  present  analysis  uses 


the  Newmark  method. 
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GOVERNING  EQUATIONS 

A  relatively  thick  rectangular-section  beam  of  thickness  h  and  length 
i  is  considered  to  be  undergoing  small -displacement  motion.  The  x  axis 
is  located  on  the  beam  midplane  and  the  z  axis  is  directed  downward  normal 
to  x. 


The  Timoshenko-beam- theory  displacement  field,  as  follows,  is  used: 
U(x,z,t)  =  u(x,t)  +  Z'Jj(x,t) 

(1) 

W(x,z,t)  *  w{x,t) 

Here,  U  and  W  are  the  displacement  components  in  the  respective  x  and  z 
directions,  u  and  w  are  the  corresponding  midplane  displacements,  t  is  time, 
and  <ii  is  the  bending  slope. 

The  axial  and  transverse  stress  resultants  N  and  Q  and  stress  couple 
M,  each  per  unit  width,  are  defined  as 


fh/2  rh/2 

(N,Q)  =  (a  ,  t  )dz  ;  M  =  Za  dz  (2) 

J-h/2  x  xz  J-h/2  x 

where  ax  and  tx2  are  respectively  the  axial  normal  stress  and  the  trans¬ 
verse  shear  stress. 

Due  to  the  bimodular  action  (different  properties  in  tension  and 
compression),  bending-stretching  coupling  is  induced  in  a  straight  beam 
of  bimodular  material  even  when  it  undergoes  small  deflections.  Thus, 
the  appropriate  constitutive  relation  is  reminiscent  of  that  for  a  beam 
unsymnetrically  laminated  of  ordinary  (not  bimodular)  materials"5: 


M 

►  = 

A  B 

r  3u/ ?x 

M 

4 

B  D 

1  3X 

^  - 

q  =  S[s|>  +  ( 3w/ 3x) ] 
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The  symbols  A,  B,  D,  and  S  denote  the  respective  stretching,  bending¬ 
stretching  coupling,  bending,  and  transverse  shear  stiffnesses,  all  per 
unit  width,  defined  by 


(A.B.D) 


rh/2 

-h/2 


Q^(l  ,z,z2)dz 


S  =  r 


th/2 

G  dz 

-h/2 


(4) 


2 

Here,  G  is  the  shear  modulus  in  the  xz  plane;  K  is  the  shear  correction 
coefficient;  Q^)  is  E^  for  a  compact-section  beam  or  E^/O  -  v  ) 

A  A  Aj  j  A 

for  a  wide  beam;  and  k=t  for  tensile-strain  regions  and  k=c  for 
compressive  regions. 

The  equations  of  motion  appropriate  for  the  subject  beam  are 
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Here,  q(x,t)  is  the  transverse  distributed  loading  and  P , I  denote  the 
respective  transverse  translatory  and  rotatory  inertia  coefficients  per 
unit  length: 

fh/2  2 

(P,I)  -  (l,z  )  dz 

-h/2 


(6) 


where  p  is  the  material  density. 

Equations  (3)  and  (5)  can  be  combined  to  obtain  the  following 
"displacement  equations  of  motion": 


_a_ 

ax 


h  Cs(f  *  *l]  ‘ p  If  - 


(7) 


As  demonstrated  in  Appendix  A,  in  bimodular  beams,  the  beam  stiff¬ 
nesses  A,  8,  and  D  depend  upon  the  neutral -surface  position  z  ,  which  is 
defined  by  setting  the  axial  normal  strain  equal  to  zero: 


Thus, 


2  o 

n  3x 


s  _  au m 

n  3x'  3x 


(8) 


TRANSFER-MATRIX  FORMULATION 

As  a  result  of  the  present  investigators’  recent  success  in  obtaining 
very  accurate  results  for  staticly  loaded,  highly  bimodular  beams  by  using 
the  transfer-matrix  method2,  it  was  decided  to  use  this  spatial-discretization 
approach  here.  Briefly,  the  beam  is  divided  into  a  number  (N$)  of  mass 
elements,  each  of  which  is  assumed  to  have  all  of  its  mass  concentrated  at 
its  mass  center,  the  location  of  which  is  called  a  station.  Consecutive 
stations  are  separated  by  massless  fields  containing  all  of  the  stiffnesses 
present  in  the  system.  Consequently,  the  complete  beam  is  discretized  to 
consist  of  two  half  fields  (di/2),  one  at  each  end  of  the  beam,  and  Ng 
stations  separated  by  Ns  -  1  full  fields  of  length  M.  The  role  of  the 
transfer  matrix,  which  will  be  presented  subsequently,  is  to  transfer  the 
generalized  displacements  (u,w,-.ji)  and  the  generalized  forces  (N,Q,M)  from 
the  left  side  of  the  field  or  station  to  the  right  side  of  the  same  field 
or  station.  For  completeness,  the  field  matrix  for  the  present  problem  is 
given  in  Appendix  B. 

Following  Chu  and  Pi 1  key9 ,  we  use  the  Newmark  acceleration  and 
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Equations 

(10)  and  (11) 

can  be  written 

in 

matrix 

notation  , 

as 
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NUMERICAL  RESULTS 

Since  no  solution,  analytical  or  numerical,  is  available  for  trans¬ 
ient  response  of  bimodular  beams,  the  present  transfer-matrix  model  could 
be  checked  only  by  comparison  with  results  obtained  for  an  ordinary- 
material  beam.  The  comparison  is  made  with  a  modal  solution  presented 
in  Ref.  9.  The  material  properties  and  dimensions  of  the  beam  are  listed 
in  Table  I  and  the  loading  and  boundary  conditions  are  those  of  Case  1 
in  Table  II.  The  tip-deflection  response  for  1 -mi  1 1 i sec  time  increments 
is  shown  in  Figure  1.  As  can  be  seen,  the  agreement  is  close,  thus 
serving  as  a  benchmark  validating  the  accuracy  of  the  present  transfer- 
matrix  model  as  applied  to  transient  analysis  of  beams. 
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The  material  and  beam  geometry  selected  for  the  bimodular  problems 
investigated  here  were  the  same  as  those  used  in  Ref.  2;  see  Table  III. 

It  is  noted  that  this  material  is  highly  bimodular  in  nature,  since 
Et/Ec  -  300. 

Transient  response  of  a  clamped-clamped  bimodular  beam  with  axial 
constraint  is  investigated  for  loadings  having  different  time  functions 
and  different  spatial  distributions.  Load  which  is  a  step  in  time  is 
investigated  in  Case  2  for  uniform  spatial  distribution  and  in  Case  3  for 
cosinusoidal  spatial  distribution.  Load  which  is  a  ramp  function  in 
time  is  investigated  in  Cases  4  and  5  for  uniform  and  cosinusoidal  dis¬ 
tribution.  Maximum-deflection  response  with  time  for  the  above  mentioned 
cases  are  illustrated  in  Figures  2  through  7.  In  all  of  these  cases, 
oscillation  keeps  continuing  without  reducing  amplitude,  since  the  model 
developed  in  this  study  does  not  incorporate  structural  damping.  However, 
the  present  model  is  still  of  value  to  designers  since  the  predicted 
results  are  on  the  conservative  side. 

CONCLUDING  REMARKS 

In  this  study,  a  transfer-matrix  model,  based  on  Timoshenko-beam 
theory,  was  developed  for  a  bimodular  beam  and  applied  to  transient 
response  problems.  Since  no  transient  response  analyses  of  this  type  of 
beam  by  other  methods  of  solution  are  presently  available,  it  was  not 
possible  to  make  comparisons.  However,  for  the  case  of  a  beam  of  ordinary 
material  subjected  to  a  triangular  pulse,  results  obtained  by  the  Dresent 
method  agreed  well  with  those  obtained  by  the  modal  method. 


minuu-a. 
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APPENDIX  A:  COMPUTATION  OF  BEAM  STIFFNESSES  FOR  A  BIMODULAR  BEAM 

The  beam  cross  section  is  assumed  to  be  compact  and  rectangular, 

(k) 

extending  from  z  =  -h/2  to  z  =  h/2.  Since  E'  '  depends  upon  the  sign  of 
the  axial  strain,  beam  stiffnesses  A,  B,  and  D  depend  upon  the  neutral- 
surface  position  zn,  thus  the  first  three  of  equations  (4)  must  be 
integrated  piecewise.  For  example,  for  concave  upward  bending 
( 32w/3x2  <  0)  the  top  portion  of  the  cross  section  (from  -h/2  to  zn)  is 
in.  compression  and  the  bottom  portion  (from  zn  to  h/2)  is  in  tension. 
Thus,  the  first  of  equations  (4)  becomes 
zn  h/2 

A  ■  Ec  dz  +  |  Et  dz  =  [zn  +  (h/2) ]EC  +  [(h/2)  -  zjE1 

-h/2 

or 

A  =  (h/2)(Et  +  Ec)  -  zn(Et-  EC) 

In  similar  fashion 

B  =  [(h2/8)  -  (zn2/2)](Et-  Ec) 

D  =  (h3/24)(Et  +  EC)  -  (zn3/3)(Et- EC) 

In  the  case  of  convex  upward  bending  (s2w/3x2  >  0),  the  results  are 
A  =  (h/2)(Et+Ec)  +  zn(Et-  Ec) 

B  -  -[(h2/8)  -  (zn2/2)](Et-EC) 

D  =  (h3/24)(Et  +  EC)  +  (zn3/3)(Et  -  EC) 
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APPENDIX  B:  FIELD  MATRIX  FOR  A  THICK  BEAM  OF  BIMODULAR  MATERIAL 
The  following  matrix  equation  was  derived  in  Ref.  2: 


2  4  6  8  10  12  14  16  Ifi  20 

t  x 10  3 ( S  ec ) 


Figure  1.  Comparison  between  transfer-matrix  and  modal  solutions 
for  a  beam  made  of  ordinary  material  and  subjected  to 
transient  load  (Case  1).  Note:  1  in.  =  2.54  cm. 


Note:  1  in.  =  2.54  cm. 


Maximum  deflection  at  (x  =  t/2)  versus  time  for 
transiently  loaded  aramid-cord  rubber  beam  (Casi 


65 


Table  I.  Material  properties  and  beam  dimensions 
for  ordinary-material  check  problem 

_ Value _ 

Property  or  Dimension  ST  units  English  units 


Material  Properties: 


Young's  modulus 
Shear  modulus 
Specific  gravity 

Length 

Depth 


2Q7  GPa  30  x  10°  psi 
79.6  GPa  11.5  x  10S  psi 

8.20  8.20 

Seam  Dimensions: 

183  cm  72.0  in. 

58.7  cm  23.1  in. 


Table  II.  Boundary  conditions  and  spatial  and  temporal 
variations  in  loading  considered 
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Table  III.  Material  properties  and  beam  dimensions 
for  aramid-cord  rubber  beam 


Value 

Property  or  Dimension 

SI  units 

English  units 

Material 

Properties : 

Longitudinal  Young's  modulus 
in  tension 

3.58  GPa 

519,300 

psi 

Longitudinal  Young's  modulus 
in  compression 

12.0  MPa 

1 ,740 

psi 

Longitudinal-thickness  shear 
modulus  (independent  of 
tension  or  compression) 

3.70  MPa 

537 

psi 

Specific  gravity 

1 .02 

1.02 

Beam  Dimensions: 

Length 

20.3  cm 

8.00 

in. 

Depth 

1 .52  cm 

0.60 

in. 
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